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ABSTRACT 


In an effort to isolate the essential physics of the meso- 
sphere and the thermosphere, a steady one-dimensional density and 
pressure model has been developed in support of related NASA 
activities, i.e., projects such as the AOTV and the Space Station. 
The model incorporates a zeroth order basic state including both 
the three-dimensional wind field and its associated shear struc- 
ture, etc. A first order wave field is also incorporated in 
period bands ranging from about one second to one day. Both 
basic state and perturbation quantities satisfy the combined 
constraints of mass, linear momentum and energy conservation on 
the midlatitude beta plane. A numerical (iterative) technique is 
used to solve for the vertical wind which is coupled to the 
density and pressure fields. The temperature structure from 1 to 
1000 km and the lower boundary conditions are specified using the 
U. S, Standard Atmosphere 1976. Vertical winds are initialized at 
the top of the Planetary Boundary Layer using Ekman pumping values 
over flat terrain. The model also allows for the generation of 
waves during the geostrophic adjustment process and incorporates 
wave nonlinearity effects. 

Preliminary results indicate that lower atmosphere wave pro- 
cesses can account for much of the observed variability of the 
density and pressure in the mesosphere and lower thermosphere. 

Basic state processes, especially vertical gradients of vertical 
winds are also influential under certain conditions. At very high 
heights (greater than 150 km) the model predicts that air density 
and pressure can deviate by as much as 20 % from Standard Atmosphere 
conditions in the presence of strong updrafts. This kind of effect 
would be likely to occur in the high latitude thermosphere where 
Joule heating is a factor leading to the development of strong 
vertical winds. 
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INTRODUCTION 


The mesosphere extends from an altitude of approximately 
50 km to approximately 90 km. The temperature in this region drops 
from about 270 K to around 180 K, and the density ranges from 
10“3 kg/m^ to 10"5 kg/m^. The mesopause separates the mesosphere 
from the thermosphere, where the temperature rises to more than 
1000 K at altitudes of hundreds of kilometers. The density of the 
thermosphere decreases to values near 10 "^ kg/m^. 

It is important for scientists and engineers working on the 
Space Shuttle, Space Station, and Aeroassisted Orbital Transfer 
Vehicle (AOTV) projects to be able to predict the air densities 
which these vehicles would encounter in the upper atmosphere, since 
large changes in the density can drastically affect their status 
while on orbit or during (Shuttle and AOTV) re-entry. In an effort 
to improve understanding of conditions in the upper atmosphere 
which can result in these changes. Dr. Douglas ReVelle has de- 
veloped a model (ReVelle, 1985) which contains some of the 
atmospheric drivers of these density variations: winds and waves. 
The computer model, called "AIR", has performed admirably at 
generating reasonable values of physical quantities. 
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OBJECTIVES 


The objectives of this work were to: 

1) transfer the program AIR from HP-75 BASIC to Applesoft BASIC 
so that 

a) editing would be easy, 

b) variables could designated with letters and numbers more 
closly resembling their names (HP-75 allows at most one 
letter and one number), 

c) graphing could be done, 

d) copies of the program could be distributed to others, 

2) run the program for a large number of different situations to 

a) determine whether or not the predictions of the model are 
reasonable, 

b) narrow choices of input parameters , 

3) modify the program to allow the input of latitude, making the 
model more "global" (calculations are now done only at 40° lati- 
tude) , 

4) rewrite the code, separating the various modes and repeated 
blocks of calculations into subroutines to make the program run 
faster, 

5) translate the program into FORTRAN, 
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THEORY 


Consider an atmosphere which is assumed to obey the following 
relationships: 

1) perfect stratification, i.e., density and pressure 
are functions of altitude only, 

2) steady state, i.e., 9/8t of any variable vanishes, 

3) inviscid flow, 

4) ideal gas law. 

The mass conservation equation is 

V-(pv) = 0 , (1) 


and the momentum equations are 


(v‘V)v + 2ftxv = -(l/p)Vp - g k > (2) 

-)■ 

where p is the density, v is the air flow velocity, is the 
angular velocity of the earth, p is the pressure, and g is the 
acceleration due to gravity. The unit vector in the vertical 
direction is £; j points north, and i points east. The wind 
velocity components are 


v x = u = zonal wind. 


v 


y 


v = meridional wind. 


v = w = vertical wind, 
z 


The theoretical model developed by ReVelle is along the lines 
of the work done by Ghosh (Ghosh, 1970) . Ghosh considered a 
steady state atmosphere with a vertical wind increasing in strength 
with height, z, and no zonal or meridional wind. Effects of 
viscosity and rotation of the earth were neglected. The solution 
of the continuity and momentum equations in this case yields a 
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density of 


p = p Q exp{-/(mg/kT + (l/T)3T/3z + (m/kT) (mg/kT 

+ (l/T)3T/3z)w 2 ) dz} , (3) 

where m is the mass of a molecule, k is Boltzmann's constant, 
and T is the temperature. The first term in the integral 
represents the gravity dependence of the density, the second term 
the temperature variation, and the third term the contribution of 
the vertical wind. 

Substitution of typical values at 100 km yields (Ghosh, 1970) 

-4 -1 

g-term ^ 1.6(10) m , 

T-term ^ 1.4(10) ** m * , 


-4 -1 

total of non-wind terms 1.7(10) m , 
-9 2 -1 

w-term 2.8(10) w m 


Therefore, if the vertical wind w = 100 m/sec. 


w-terms 
non w-terms 


^ 10 %. 


ReVelle considers a rotating earth with non-zero zonal wind 
and meridional winds and non-zero spatial derivatives of the 
winds. The continuity equation becomes 


p3u/3x + p3v/3y + p3w/3z + w3p/3z = 0 , (4) 

or pD^ + p3w/3z + w3p/3z = 0 , (5) 

where D^, the horizontal divergence, is assumed to be constant 
for all z. 
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The atmosphere is broken up into layers of thickness Az (^.1 km), 
and the conservation equation and momentum equations are numeri- 
cally integrated from a lower boundary of altitude z = 1 km. The 
initial value of the vertical wind is taken from Ekman pumping 
theory in which 


w 

o 


^ C z h PBl/ 6 


( 6 ) 


where C z is the vorticity (in sec ) of some rotating body of air 
at the surface and hpg L is the height of the planetary boundary 
layer, 1 km. 

Note that if w is assumed a constant, w 0 , over some layer of 
thickness Az, then equation (5) can be written as 


p w = P o w o exp{-CD h /w o )Az} . (7) 

The case in which = 0 and the case in which w is very large 
should both act similarly with pw 'v* constant. Also note that if 
w becomes negative (a downdraft), the exponent becomes positive, 
causing the density to increase with height. This sign in the 
exponent prevents a downdraft from becoming an updraft again; 
therefore, the solutions which were taken to be acceptable in 
this study were the ones in which the vertical wind was always 
upward . 

It would be desirable to also utilize the conservation of 
energy equation in the model. However, due to the complexity of 
the equation and the uncertainty in values of some variables, it 
was decided to stick with only the mechanical relationships. A 
rough calculation does reveal that the winds that would be 
expected from an energy point of view are compatible with those 
being calculated by the model. The heating due to vertical 
lifting is calculated in the program (variable name Q7) , and 
the numbers obtained are consistent with calculations which show 
heating and cooling rates of 1 - 10 K/day in the lower 100 km 
of atmosphere and up to 1000 f s K/day in the thermosphere. 

The details of using the program are given in the next 
section. The rest of this section is devoted to a discussion of 
a few items essential to the understanding of the operation of 
the program. 

There are a number of program input requirements; the first 
is that of mode. The program operates in one of four modes: 
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1) barotropic atmosphere with vertical winds (the only mode 
run in this study) 

2) sub/super geostrophic with either constant or height 
variable vertical winds) 

3) sub/super geostrophic with waves (includes nonlinear 
saturation limit) 

4) hydrostatic with waves (includes nonlinear saturation 
limit) . 


A response of "yes" to the input concerning heating constraint 
forces the vertical wind to become constant when the heating due 
to vertical lifting becomes greater than 100 K/day. A response 
of "no" allows the vertical wind and its associated heating to 
increase without bound. 

The value of vorticity input is used to calculate the initial 
vertical wind in equation (6) . The number most often used in 
this study is 0 .606433553(10) sec"l; it generates an initial 
vertical wind of 1.00687771 mm/sec. This vorticity was found by 
ReVelle prior to this study to yield results (winds, densities, 
etc.) which are realistic and well behaved. 

The DEL:H DOT V and DELH.V terms in the program refer to the 
horizontal divergence, = 9u/9x + 9v/9y, in equation (5) . A 
response of "yes" to the "NOMINAL DELH.V CASE?" input causes the 
value of Dfo to be set at 1.21915 (10) sec - *, which was also 
previously determined by ReVelle to yield reasonable results. 

An input of "no" prompts a request for the desired value of Dft. 
Whatever the value of D^ used in the program, it is held constant 
everywhere . 

The zonal wind can be either constant with altitude or 
variable according to the season. A constant zonal wind is ob- 
tained by answering "no" to the "VARIABLE CORIOLIS TERM?" input 
and then typing in the desired wind. Otherwise, the zonal wind 
is variable according to whether winter or summer is chosen. In 
this case the wind is assumed geostrophic such that 


u = - (R/fM) 9T/9y , (8) 

where f is the Coriolis term (2£lsin0) , R is the gas constant, M 
is molecular mass per mole, and 9T/3y is the horizontal temperature 
gradient. The zonal wind may in addition be made subgeostrophic. 
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geostrophic, or supergeostrophic via the u geostrophic factor 
(<1, 1, >1, respectively) which simply multiplies the value 
calculated in equation (8) . 


OPERATING INSTRUCTIONS 

Step-by-step instructions on how to use the Apple 11+ or 
lie (or III with II emulator) to run AIR and its derivatives are 
given in this section. 

1) Insert the disk into the disk drive and close it. 

2) Turn on the computer (switch is on your left, at rear of 

computer) and the monitor. 

3) The menu of the disk contents will appear on the screen. 

4) Type RUN AIR ( or RUN AIR PRINT or RUN AIR LARGE SCALE) 

and press the "return" key. 

5) Respond to each input request by typing the number and 
pressing "return". 

6) When the run is complete, remove the disk and turn off 
the computer and monitor. 

Notes on AIR: 

* Input 

MODE (1, 2, 3, 4) = 

(Mode 1 was used for all runs.) 

HEATING CONSTRAINT: 1(Y), 0(N) = 

(0 used for all runs.) 

PBL Z VORTICITY: MEAN (l*E-5/SEC) 

(0.606433553 used for most runs.) 

NOMINAL DELH.V CASE? (1/YES, O/NO) 

(YES chooses 1.21915.) 

VARIABLE CORIOLIS TERM? (1/YES, O/NO) 

(Variable or constant zonal wind.) 
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ZONAL WIND = 

(Appears only if Variable Coriolis Term? = NO.) 

WINTER/SUMMER CASE? (1/WIN, 0/SUM) 

(Appears only if Variable Coriolis Term? = YES.) 

U GEOSTROPHIC FACTOR = 

(1 used for all runs.) 

DEL:H DOT V(l*E-7/SEC) = 

Horizontal divergence in units of (10)" sec" ; usually 
in the range of 1.2. (Appears only if Nominal DEL.V 
Case? = NO.) 


* Output 

Any number of variables can be listed on the screen. 
Which variables are to be listed is set in line number 18700 
(and 18710 if the list is very long) . The altitude, L8, is 
already part of the PRINT statement. If, for example, the 
variable Y4 is also needed, then this sequence of steps 
should be followed: 

1) Instead of typing RUN AIR, type LOAD AIR. 

2) Type LIST 18700-18800 (press "return") to see what the 
PRINT statement looks like. 

3) Retype line 18700. In this example it would be 
18700 PRINT :PRINT "L8 = ";LR : PRINT "Y4 = ";Y4 

4) Type RUN and press "return". 

The program may be stopped at any time by pressing the CTRL 
("control") and C keys simultaneously. The current value of 
a variable, for example R6, is listed by typing ?R6 (press 
"return") . The program operation is continued with CONT. 

* Some useful variables and their code names : 

LR - altitude (rounded to the nearest 0.01 km) 

T1 - vertical temperature gradient 
T2 - temperature 
M - density 

M9 - standard atmosphere density 
PI - pressure 

P9 - standard atmosphere pressure 
Y4 - % difference between M and M9 
Y9 - % difference between PI and P9 
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R6 - Reynolds number 
G1 - density scale height 
S2 - pressure scale height 
W - vertical wind 
MO - Mach number 

Q7 - heating due to vertical lifting 

U2 - zonal wind 

HO - horizontal divergence 


Notes on AIR LARGE SCALE: 

This program was created to eliminate the need to modify 
AIR every time a run up to 300 km was desired. The main 
difference between AIR LARGE SCALE and AIR is in the defini- 
tion of the x- and y-coordinates on the graph. The altitude 
scale goes to 300 km (as opposed to 150 km in AIR) and the 
wind scale goes to 14 mm/sec (2.8 mm/sec in AIR). 

Other variables may be graphed by redefining the array 
variable W(Z) used in the graphics. W(Z) is dimensioned in 
line 50 and defined in line 18630. 


Notes on AIR PRINT : 

This program was created to give a paper printout of the 
variables being listed in line 18700. The printer is accessed 
through slot #1 (PR#1 in line 3310), and it must be turned on 
before the program is run. The program is set to run until 
the altitude reaches 500 km, but it can be stopped at any 
point with CTRL-C. Type PR#0 after this to return to use of 
the monitor. 


RESULTS 

This section comprises a summary of the results of the runs 
of AIR and their significance. All runs were made with the 
vorticity = 0. 606433553(10) “5 sec - l unless otherwise stated. 

1) u = 10 m/sec, nominal 

Figure 1 shows a typical behavior of the vertical wind. A 
later run of this case to higher heights shows that the vertical 
wind increases to 10 m/sec at 155 km, where the density is 0.2 % 
greater than the 1976 Standard Atmosphere (National Oceanic and 
Atmospheric Administration, 1976) density, to 100 m/sec at 214 km. 
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where the density is 1.3 % less than standard, to almost 400 m/sec 
at 257 km, where the density drops to 16.6 % less than standard. 

It is noted here that the shape of the curve in Figure 1 is 
very similar to the temperature profile. Figure 2. Because of this 
similarity, it was decided to attempt to find values of and u 
which would cause the w vs. z graph to follow the temperature 
curve. The results of this search are given later in this section. 

2) u = 0, nominal 

Figure 3 shows how the vertical wind profile changes with the 
zonal wind reduced to zero. A negative (toward the west) zonal 
wind drives the vertical wind to large values at even lower 
heights. See Figure 4. 

-7 -1 

3) u = 10 m/sec, = 1.2(10) sec and 

u = 10 m/sec, = 1.24(10)"^ sec 

These graphs demonstrate how sensitive the w vs. z graphs are 
to the value of the horizontal divergence. A value of only 
slightly smaller than 1.21915(10) sec - * causes the wind to in- 
crease very rapidly (Figure 5) . In fact a zero divergence yields 
a vertical wind of 200 m/sec at only 70 km. A value of slightly 
larger than nominal drives the wind negative very quickly. See 
Figure 6. 

A comparison of the u = 0 and = 0 case with Ghosh's re- 
sults might be in order, although they will still differ somewhat 
due to Ghosh's assumption of a non-rotating earth. A run with 
u = 0 and = 0 yields a vertical wind of 100 m/sec at 83 km 
(y 100 km), which produces a density 9 % (y 10 %) smaller than 
Standard Atmosphere density, which includes the gravity and 
temperature variations . 

4) u = 11.723013 m/sec, nominal 

Figure 7 shows that this value of a constant zonal wind 
constrains the wind to follow the temperature profile for a good 
distance and limits it to a few cm/sec near 300 km. 

Figure 8 demonstrates how the winds vary as a result of a 
change in only the last digit of the zonal wind. The vertical 
wind is just beginning to decrease at 300 km and will eventually 
go negative. 

5) Winter and summer zonal winds 

Figures 9 and 10 show the zonal (geostrophic) wind profile 
for winter and summer, respectively. 
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w(mm/sec) 

Figure 1. Vertical wind vs. altitude for u = 10 m/sec 
and D h = 1 . 21915 (10) " 7 sec" 1 . 



T(K) 

Figure 2. Temperature profile. 
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Figure 3. Vertical wind vs. altitude for u = 0 m/sec 
and D h = 1 .21915 (10) " 7 sec -1 .' 



Figure 4. Vertical wind vs. altitude for u = -10 m/sec 
and D h = 1.21915(10)" 7 sec' 1 . 

XIII-12 



300 



Figure 7. Vertical wind vs. altitude for u = 11.723013 m/sec 
and D h = 1 .21915 CIO) " 7 sec -1 . 



Figure 8. Vertical wind vs. altitude for u = 11.723014 m/sec 
and D h = 1.21915(10)" 7 sec' 1 . 
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Figure 10. Summer zonal wind profile. 
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6) u = summer profile, = 1.21917253(10) sec 

This value of the horizontal divergence was found to be the 
one which constrained the vertical wind to follow the temperature 
profile as closely as possible in the case of summer zonal winds. 

The vertical wind (see Figure 11) is only about 1 cm/sec at 170 km 
and rises to 10.4 m/sec at 500 km, where the density is 3.2 % 
above the Standard Atmosphere density. 

7) u = winter profile, = 1 ,21892828(10) " 7 se c 1 

Similarly, this constraining value of was found for the 
winter case. See Figure 12. The vertical wind increases to 
1 cm/sec at 155 km, where the density is 1.3 % greater than 
standard, to 21.5 m/sec at 500 km, where the density is 4.9 % 
greater than standard. 

8) u = 10.105 m/sec, D h = 1.21917253(10)“ 7 sec" 1 

It was noted early in the study that for a nominal the 
w vs. z graphs for u = 10 m/sec and u = summer profile were almost 
identical and that the u = 25 m/sec and u = winter profile were 
almost identical. The ratio of these two constant zonal winds is 
2.5, which is roughly equal to the ratio of the altitude -averaged 
zonal winds for winter and summer. Therefore, once the constraining 
values of the horizontal divergence were found for winter and 
summer, the corresponding constant zonal winds for these divergences 
were sought. Figure 13 shows how the vertical wind follows the 
temperature profile for u = 10.105 m/sec. 

9) u = 27.6 m/sec, D h = 1.21892828(10)" 7 sec" 1 

The constrained vertical wind for the case of a constant 
zonal wind and the winter divergence is graphed in Figure 14. 

Again, note that the ratio of these two constant zonal winds in 
the winter and summer cases is 2.73 (y 2.5). Time did not permit 
the determination of the zonal winds to more decimal places. 

10) Note on heating 

As discussed earlier, if the theory is to be reasonable, the 
heating due to the vertical winds must be within the limits men- 
tioned. A run of the constrained summer case produced vertical 
wind heating which a) followed the temperature profile, b) was 
less than or approximately equal to 1 K/day up to 108 km, c) hit 
100 K/day at 279 km, and d) rose to 3700 K/day at 500 km. 

i. 

The heating due to vertical lifting in the constrained winter 
case behaved much like the summer case except that the heating hit 
100 K/day at 242 km and 7700 K/day at 500 km. 
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Figure 13. 


Vertical 
and Dh = 


wind vs. altitude for u = 10.105 m/sec 
1.21917253(10)'” 


see"* 





Figure 14. Vertical wind vs. altitude for u = 27.6 m/sec 
and = 1 .21892828(10) " 7 sec -1 . 


< 

I 


XIII-18 


11) Note on mass flux 

A run was made to test the earlier assertion that as w gets 
very big, pw, the mass flux in kg*m“^*sec“l, should become 
constant. Figures 15 and 16 give the wind vs. altitude and mass 
flux vs. altitude graphs, respectively. The vertical wind is 
435 m/sec at 263 km, and the density is 21 % smaller than 
standard. The mass flux reaches a constant value of 1.5(10)”® 
kg*m~2*sec“l a t around 120 km. 

At this rate of mass flux the earth 1 s atmosphere would be 
lost in 


5 . 3(10) 18 kg 

t -v * 5.2(10) sec 

{1.5(10)"^ kg/m^*sec}{4tr(6. 7(10)^m)^} 

or about 20,000 years. Even if it is only the upper 1 % that is 
lost, the time is reduced to 200 years. These numbers are not 
meant to show that the atmosphere is actually dissipating at 
this rate, rather that the quantities produced by AIR are reason- 
able. That is, the time is not as short as 2 seconds, nor is it 
as long as 2 billion years. 

12) Note on pressure variations 

All runs yielded a pressure variation (variable Y9) very 
close to the density variation (variable Y4) . So, when the den- 
sity was a certain percent above (below) the standard, the pressure 
was also that percent above (below) the standard. 


XIII-19 



-2 


log(wCm/sec)) 


logCpwCkg^nT^ •sec - -*-)) 


Figure 15. 
tude for z, 
10 ra/ sec, 


Vertical wind vs. alti- 
z = 0.57(10) " 5 sec" 1 , u « 
Dh = 1 ,21915 (10) " 7 sec -1 . 


Figure 16. Mass flux vs. alti- 
tude for c z = 0.57(10)~S sec -1 , u = 
10 m/sec, D h = 1 . 21915 (10) -7 sec -1 . 
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CONCLUSIONS AND RECOMMENDATIONS 


Objectives 1) and 2) were achieved with the possible exception 
of lb); some variables have been renamed, but the changes have 
not been made in the code. Due to the fact that the code was still 
in a state of flux throughout most of the summer, time was short 
and objectives 3) and 4) were not achieved. They will remain as 
long term objectives. Objective 5) was never really attempted 
because a) the Apple microcomputer is so convenient in terms of 
accessibility and editing and graphing capabilities, b) time on 
the bigger computers is very limited, c) a decrease in run time 
by a factor of four on the Hewlett -Pachard 1000F is, although 
desirable, not significant, and d) time, as mentioned above, was 
limited. Perhaps, once the present code is brought to its final 
form, the translation to FORTRAN can take place. 

As a result of this study, the following recommendations are 
put forth: 

1) continue work on the code to achieve objectives 3) and 4), 

2) extend the work done on objective 2) to the other modes 
in the program, 

3) study the possibility of allowing for a variable 
horizontal divergence, rather than keeping it constant 
for all altitudes. 
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